function P_kp1=Nonlin(Sigma, dSigma, dTau, N_Nonlin, N_t, N_r, P_k)

Dummy=N_Nonlin*log(1 + dSigma/(1 + Sigma));
P_kp1=zeros(size(P_k));

for j=1:N_r - 1
  for i=2:N_t - 1
    if P_k(i, j)>=0
      P_kp1(i, j)=P_k(i, j)/(1 - Dummy*(P_k(i + 1, j) - P_k(i, j))/dTau);
    else
      P_kp1(i, j)=P_k(i, j)/(1 - Dummy*(P_k(i, j) - P_k(i - 1, j))/dTau);
    end  
  end  
end
% note: Note that 'N_r' runs from '1' to 'N_r - 1', which is required for
% the correct handling of the boundary conditions. Including 'N_r' will
% however, not result in significant errors in most cases.
% note: Strang enough the double for-loop implementation is also the
% fastest one. A matrix implementation is more than 2 times slower.